clc
clear
rho1=1000;rho2=600;rho3=200;
h1=5000;h2=5000;
miu=4*pi*10^(-7);
omega=logspace(-4,4,100);
k1=sqrt(-i*omega*miu/rho1);
k2=sqrt(-i*omega*miu/rho2);
k3=sqrt(-i*omega*miu/rho3);
R2=tanh((h2.*k2+atanh(k2./k3)));
R1=tanh(h1*k1+atanh((k1./k2).*R2));
y=rho1*(R1).^2;
figure(1)
semilogx(omega,abs(y));
title('Q型振幅曲线')
xlabel('频率')
ylabel('视电阻率')
ang=2*angle(R1)*180;
figure(2)
semilogx(omega,ang);
title('Q型相位曲线')
xlabel('频率')
ylabel('相位')

% 第三层电阻率为无限大时代码：
clc
clear
rho1=900;rho2=1/9*rho1;rho3=inf;
a=1/3;
h1=10000;
h2=h1*a;      %第一层厚度和第二层厚度
miu=4*pi*10^(-7);    %真空介电常数
omega=logspace(-9,6,100);
k1=sqrt(-i*omega*miu/rho1);
k2=sqrt(-i*omega*miu/rho2);
k3=sqrt(-i*omega*miu/rho3);
y1=h1.*k1;
y2=h2.*k2;
R2=tanh(h2.*k2+atanh(k2./k3));
y=rho1*tanh(y1+atanh((k2./k1).*R2)).^2;
loglog(1./omega,abs(y)/rho1,'linewidth',2);
hold on
 %%
a=1;
h2=h1*a;
miu=4*pi*10^(-7);
 
k1=sqrt(-i*omega*miu/rho1);k2=sqrt(-i*omega*miu/rho2);k3=sqrt(-i*omega*miu/rho3);
y1=h1.*k1;
y2=h2.*k2;
R2=tanh(h2.*k2+atanh(k2./k3));
y=rho1*tanh(y1+atanh((k1./k2).*R2)).^2;
loglog(1./omega,abs(y)/rho1,'linewidth',2);
 %%
a=5;
h2=h1*a;
y1=h1.*k1;
y2=h2.*k2;
R2=tanh(h2.*k2+atanh(k2./k3));
y=rho1*tanh(y1+atanh((k1./k2).*R2)).^2;
loglog(1./omega,abs(y)/rho1,'linewidth',2);
 %%
a=9;
h2=h1*a;
y1=h1.*k1;
y2=h2.*k2;
R2=tanh(h2.*k2+atanh(k2./k3));
y=rho1*tanh(y1+atanh((k1./k2).*R2)).^2;
loglog(1./omega,abs(y)/rho1,'linewidth',2);
%%
a=24;
h2=h1*a;
y1=h1.*k1;
y2=h2.*k2;
R2=tanh(h2.*k2+atanh(k2./k3));
y=rho1*tanh(y1+atanh((k1./k2).*R2)).^2;
loglog(1./omega,abs(y)/rho1,'linewidth',2);
%%
legend('1/3','1','5','9','24')
title('底层电阻率区域无穷');
ylabel('\rho_\omega\\\rho_1');
xlabel('1\\\omega')
